PATENT 



vo Attorney Docket No. CANNING. 001 A 

^=jSJ Date: September 29, 2000 

SB c Pa § el 

^5 2 ASSISTANT COMMISSIONER FOR PATENTS o 
WASHINGTON, D.C. 20231 

CM 

ATTENTION: BOX PATENT APPLICATION 

Sir: ^ 
Transmitted herewith for filing is the patent application of u*" 3 
Inventor(s): Francis X. Canning 

For: COMPRESSION AND COMPRESSED INVERSION OF INTERACTION DATA 

Enclosed are: 

(X) Specification in thirty-six (36) pages. 

(X) Thirteen (13) sheets of drawings. 

™ (X) A verified statement to establish small entity status under 37 CFR 1 .9 and 37 CFR 1 .27. 

(X) Initial signed Declaration and Power of Attorney by inventor(s). 

%j (X) Return prepaid postcard. 



CLAIMS AS FILED 



FOR NUMBER 

FILED 


NUMBER 
EXTRA 


RATE 


FEE 


Basic Fee 




$345 


$345 


Total Claims 33 - 20 = 


13 x 


$9 


$117 


Independent Claims 6 - 3 = 


3 x 


$39 


$117 


If application contains any multiple dependent claims(s), then add 


$130 


$0 


TOTAL FILING 
FEE 


$579.00 







(X) A check in the amount of $579.00 to cover the filing fee is enclosed. 



(X) The Commissioner is hereby authorized to charge any additional fees which may be required, now or in the 
future, or credit any overpayment to Account No. 1 1-1410. A duplicate copy of this sheet is enclosed. 



(X) Please use Customer No. 20,995 for the correspondence address. 




Lee W. Henderson Ph.D. 
Registration No. 41,830 
Attorney of Record 

H:\DOCS\LWH\LWH-4827.DOC:pp/092900 



KNOBBE, MARTENS, OLSON & BEAR, LLP 
620 NEWPORT CENTER DR 16TH FLOOR NEWPORT BEACH, CA 92660 

C949) 760-0404 FAX 760-9502 



A 



INTELLECTUAL PROPERTY LAW 

KNOBBE, MARTENS, OLSON 8c BEAR 



LOUIS J. KNOBBE 
DON W. MARTENS* 
GORDON H. OLSON* 
JAMES B. BEAR 
0 ARRELL L. OLSON* 
WILLIAM B. BUNKER 
WILLIAM H. N1EMAN 
ARTHUR S. ROSE 
JAMES F. LESNIAK 
NED A. ISRAELSEN 
DREW S. HAMILTON 
JERRY T. SE WELL 
JOHN B. SGANGA, JR. 
EDWARD A. SCHLATTER 
GERARD VON HOFFMANN 
JOSEPH R. RE 
CATHERINE J. HOLLAND 
JOHN M. CAR SON 
KAREN VOGEL WEIL 
ANDREW H. SIMPSON 
JEFFREY L. VAN HOOSEAR 
DANIEL E. ALTMAN 
MARGUERITE L..GUNN 
STEPHEN C. JENSEN 
VITO A. CANUSO III 
WILLIAM H. SHREVE 
LYNDA J. ZADRA-SYMES + 
STEVEN J. NATAUPSKY 
PAUL A. STEWART 
JOSEPH F. JENNINGS 
CRAIG S. SUMMERS 
ANNEMARIE KAISER 
BRENTON R. BABCOCK 

— u. 



THOMAS F SMEGAL, JR 
MICHAEL H TRENHOLM 
DIANE M REED 
JONATHAN A BARNEY 
RONALD J SCHOENBAUM 
JOHN R KING 
FREDERICK S B E R R ETTA 
NANCY WAYS VENSKO 
JOHN P GIE2EN TANNER 
ADEEL S AKHTAR 
GINGER R DREGER 
THOMAS R ARNO 
DAVID N WEISS 
DANIEL HART, PH.D 
DOUGLAS G MUEHLHAUSER 
LORI LEE YAMATO 
MICHAEL K FRIEDLAND 
STEPHEN M LOBBIN 
ST AC E Y R HALPERN 
DALE C HUNT, PH D 
LEE W HENDERSON, PH D 
DEBORAH S SHEPHERD 
RICHARD E CAMPBELL 
MARK M ABUMERI 
JON W GURKA 
ERIC M. NELSON 
MARK R BENEDICT, PH D 
PAUL N CONOVER 
ROBERT J ROB Y 
SABING H LEE 
KAROLINE A, OELANEY 
JOHN W HOLCOMB 
JAMES J MULLEN, III, PH D 



A LIMITED LIABILITY PARTNERSHIP INCLUDING 
PROFESSIONAL CORPORATIONS 

PATENT, TRADEMARK AND COPYRIGHT CAUSES 

620 NEWPORT CENTER DRIVE 
SIXTEENTH FLOOR 
NEWPORT BEACH, CALIFORNIA 92660-8016 
C94-SO 760-0404 
FAX 760-9502 
INTERNET. WWW KMOB COM 



JOSEPH S CIANFR ANI 
JOSEPH M REISMAN, PH.D 
WILLIAM R. ZIMMERMAN 
GLEN L NUTTALL 
ERIC S FURMAN, PH D 
TIRZAH ABE LOWE 
GEOFFREY Y MDA 
ALEXANDER S FRANCO 
SANJIVPAL S. GILL 
SUSAN M MOSS 
JAMES W HILL, M D 
ROSE M THIESSEN, PH D 
MICHAEL L FULLER 
MICHAEL A GUILIANA 
MARK J KERTZ 
RABINDER N NARULA 
BRUCE S ITCHKAWITZ, PH D 
PETER M MIDGLEY 
THOMAS S MCCLENAHAN 
MICHAEL S OK AMO TO 
JOHN M GROVER 
MALLARY K DE MERLIER 
! R F AN A LATEEF 
AMY C CHRISTENSEN 
SHARON S NG 
MARK J GALLAGHER, PH D 
DAVID G JANKOWSKI, PH D 
BRIAN C HORNE 
PAYSON J LEMEILLEUR 
WILLIAM G BERRY 
DIANA W PRINCE 



*> Assistant Commissioner for Patents 
c Washington, D.C. 20231 



i TO 



CERTIFICATE OF MAILING BY "EXPRESS MAIL" 



Attorney Docket No. 

Applicant(s) 

For 



OF COUNSEL 



JERRY R SEILER 
PAUL C STEINHARDT 



JAPANESE PATENT ATT Y 
K ATSUHIRO AR A I** 



EUROPEAN PATENT ATTY 
MARTIN HELLEBRANDT 



KOREAN PATENT ATTY 
MINCHEOL KIM 

SCIENTISTS & ENGINEERS 
(NON-LAWYERS) 

RAIMOND J SALENIEKS** 
DANIEL E JOHNSON, PH D ** 
JEFFERY KOEPKE, PH D ** 
KHURRAM RAHMAN, PH D 
JENNIFER A HAYNES, PH D. 
BRENDAN P 0 NEILL, PH D 
THOMAS Y MAG AT A 
LINDA H LIU 

YASHWANT VAISHNAV, PH D 

MEG UM I T ANA K A 

CHE S CHERESKIN, PH D ** 

ERIK W ARCHBOLD 

PHILIP C. HARTSTE1N 

JULIE A HOPPER 

CHRIS S CASTLE 

JAMES W AUSLEY 

R P CARON, PH D 

JENNIFER HAYES 

KIRK E PASTORIAN, PH D 

CHARLES T. RIOGELY 

KEITH R MCCOLLUM 

LANG J MCHARDY 

* A PROFESSIONAL CORPORATION 
t ALSO BARRISTER AT LAW <U K ) 
**US PATENT AGENT 



CANNING.001A 
Francis X. Canning 

COMPRESSION AND COMPRESSED 
INVERSION OF INTERACTION DATA 




Attorney 



Lee W. Henderson Ph.D. 



"Express Mail" 
Mailing Label No. 

Date of Deposit 



EL559446715US 
September 29, 2000 



I hereby certify that the accompanying 



Transmittal in Duplicate; Specification in 36 pages; 13 sheets of drawings; 
SIGNED Declaration and Power of Attorney in 1 pages; Small Entity Statement; 
Check for Filing Fee; Return Prepaid Postcard 



are being deposited with the United States Postal Service "Express Mail Post Office to 
Addressee" service under 37 CFR 1.10 on the date indicated above and are addressed to the 
Assistant Commissioner for Patents, Washington, D.C. 20231. 



H:\DOCS\LWH\LWH-4828.DOC :pp/092900 




201 CALIFORNIA STREET 550 WEST C STREET 3301 UNIVERSITY AVENUE 1900 AVENUE OF THE STARS 

SUITE 1150 SUITE 1200 SUITE 710 SUITE 1 + 25 

SAN FRANCISCO, CALIFORNIA 9 + 111 SAN DIEGO, CALIFORMIA 92101 RIVERSIDE, CALIFORNIA 92501 LOS ANGELES, CALIFORNIA 90057 

(415) 954-4114 (619) 235-B55D (909) 781-9231 (310) 551-3450 

FAX (415) 954-4111 FAX (619) 235-0176 FAX <909) 781-4507 FAX (310) 551-345B 



Applicant or Patentee; Francis X. Canning Attorney's Docket No.: CANNING.001A 

Application or Patent No.: Unknown Page 1 

Filed or Issued: Herewith 

For: COMPRESSION AND COMPRESSED INVERSION OF INTERACTION DATA 



VERIFIED STATEMENT (DECLARATION) CLAIMING SMALL-ENTITY STATUS 

I, the undersigned, do hereby declare that: 

[X] I am an independent inventor as defined in 37 CFR 1.9(c) for purposes of paying reduced fees to the Patent 
and Trademark Office with regard to the invention described in the patent or application identified above. 

The individual, concern or organization identified above has not assigned, granted, conveyed or licensed, and is 
under no obligation under contract or law to assign, grant, convey or license, any rights in the invention to any 
person who would not qualify as an independent inventor under 37 CFR 1.9(c) if that person had made the 
invention, or to any concern which would not qualify as a small business concern under 37 CFR 1.9(d) or a 
nonprofit organization under 37 CFR 1.9(e). 

If the rights held by the above-identified individual, concern or organization are not exclusive, each individual, 
concern or organization having rights in the invention are identified below. Each such individual, concern or 
organization must file separate verified statements averring to their status as small entities. 

*NOTE: Separate verified statements are required from each named person, concern or organization having 
rights to the invention averring to their status as small entities. (37 CFR 1.27). 

FULL NAME: Francis X. Canning 

ADDRESS: 3557 Lesser Drive, Newbury Park, CA 91320 

[X] INDIVIDUAL [ ] SMALL BUSINESS CONCERN [ ] NONPROFIT ORGANIZATION 

I acknowledge the duty to file, in this application or patent, notification of any change in status resulting in loss of 
entitlement to small-entity status prior to paying, or at the time of paying, the earliest of the issue fee or any 
maintenance fee due after the date on which status as a small entity is no longer appropriate. (37 CFR 1.28(b)). 

I hereby declare that all statements made herein of my own knowledge are true and that all statements made on 
information and belief are believed to be true; and further that these statements were made with the knowledge that 
willful false statements and the like so made are punishable by fine or imprisonment, or both, under section 1001 of 
Title 18 of the United States Code, and that such willful false statements may jeopardize the validity of the 
application, any patent issuing thereon, or any patent to which this verified statement is directed. 



NAME OF PERSON SIGNING: Francis X. Canning 

ADDRESS OF PERSON SIGNING: 3557 Lesser Drive, Newbury Park, CA 91320 
SIGNATURE: ^\^rS^P^ ^^fx DATE: <^^^fc^^ 3 % 



H:\DOCS\LWH\LWH-4680.DOC :pp/082800 




CANNING.001A 



PATENT 



COMPRESSION AND COMPRESSED INVERSION OF INTERACTION DATA 

Related Applications 

The present application claims priority under 35 U.S.C. 119(e) from U.S. Provisional 
Application No. 60/175454, filed on January 10, 2000, titled "COMPRESSION AND 
COMPRESSED INVERSION OF INTERACTION DATA" and from U.S. Provisional 
Application No. 60/201149, filed on May 2, 2000, titled "COMPRESSION AND 
COMPRESSED INVERSION OF INTERACTION DATA," the contents of which are hereby 
incorporated by reference in their entirety. 

Background of the Invention 

Field of the Invention 

The invention relates to methods for compressing the stored data, and methods for 
manipulating the compressed data, in numerical solutions involving numerous mutual 
interactions, especially when the nature of these interactions approaches an asymptotic form 
for large distances, such as, for example, antenna problems solved using the method of 
moments. 

Description of the Related Art 

Many numerical techniques are based on a "divide and conquer" strategy wherein a 
complex structure or a complex problem is broken up into a number of smaller, more easily 
solved problems. Such strategies are particularly useful for solving integral equation problems 
involving radiation, heat transfer, scattering, mechanical stress, vibration, and the like. In a 
typical solution, a larger structure is broken up into a number of smaller structures, called 
elements, and the coupling or interaction between each element and every other element is 
calculated. For example, if a structure is broken up into 16 elements, then the inter-element 
mutual interaction (or coupling) between each element and every other element can be 
expressed as a 16 by 16 interaction matrix. 



As computers become more powerful, such element-based numerical techniques are 
becoming increasingly important. However, when it is necessary to simultaneously keep track 
of many, or all, mutual interactions, the number of such interactions grows very quickly. The 
size of the interaction matrix often becomes so large that data compression schemes are 
5 desirable or even essential. Also, the number of computer operations necessary to process the 
data stored in the interaction matrix can become excessive. The speed of the compression 
scheme is also important, especially if the data in the interaction matrix has to be 
decompressed before it can be used. 

Typically, especially with radiation-type problems involving sound, vibration, stress, 
10 temperature, electromagnetic radiation, and the like, elements that are physically close to one 
another produce strong interactions. Elements that are relatively far apart (usually where 
O distance is expressed in terms of a size, wavelength, or other similar metric) will usually 
pi couple less strongly. For example, when describing the sound emanating from a loudspeaker, 
^ the sound will change in character relatively quickly in the vicinity of that speaker. If a person 
15 standing very near the speaker moves one foot closer, the sound may get noticeably louder. 
\| However, if that person is sitting at the other end of a room, and moves one foot closer, then 
ri the change in volume of the sound will be relatively small. This is an example of a general 
J! property of many physical systems. Often, in describing the interaction of two nearby objects, 
yy relatively more detail is needed for an accurate description, while relatively less detail is 
p 20 needed when the two objects are further apart. 

As another example, consider a speaker producing sound inside a room. To determine 
the sound intensity throughout that room, one can calculate the movement (vibration) of the 
walls and objects in the room. Typically such calculation will involve choosing a large 
number of evenly spaced locations in the room, and determining how each location vibrates. 
25 The vibration at any one location will be a source of sound, which will typically react with 
every other location in the room. The number of such interactions would be very large and 
the associated storage needed to describe such interactions can become prohibitively large. 
Moreover, the computational effort needed to solve the matrix of interactions can become 
prohibitive. 
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Summary of the Invention 
The present invention solves these and other problems by providing a compression 
scheme for interaction data and an efficient method for processing the compressed data 
without the need to first decompress the data. In other words, the data can be numerically 
5 manipulated in its compressed state. 

Given a first region containing sources relatively near to each other, and a second 
region containing sources relatively near to each other, but removed from the first region; one 
embodiment provides a simplified description of the possible interactions between these two 
regions. That is, the first region can contain a relatively large number of sources and a 

10 relatively large amount of data to describe mutual interactions between sources within the first 
region. In one embodiment, a reduced amount of information about the sources in the first 
region is sufficient to describe how the first region interacts with the second region. One 
embodiment includes a way to find these reduced interactions with relatively less 
computational effort than in the prior art. 

15 For example, one embodiment includes a first region of sources in one part of a 

problem space, and a second region of sources in a portion of the problem space that is 
removed from the first region. Original sources in the first region are modeled as composite 
sources (with relatively fewer composite sources than original sources). In one embodiment, 
the composite sources are described by linear combinations of the original sources. The 

20 composite sources are reacted with composite testers to compute interactions between the 
composite sources and composite testers in the two regions. The use of composite sources 
and composite testers allows reactions in the room (between regions that are removed from 
each other) to be described using fewer matrix elements than if the reactions were described 
using the original sources and testers. While an interaction matrix based on the original 

25 sources and testers is typically not a sparse matrix, the interaction matrix based on the 
composite sources and testers is typically a sparse matrix having a block structure. 

One embodiment is compatible with computer programs that store large arrays of 
mutual interaction data. This is useful since it can be readily used in connection with existing 
computer programs. In one embodiment, the reduced features found for a first interaction 
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group are sufficient to calculate interactions with a second interaction group or with several 
interaction groups. In one embodiment, the reduced features for the first group are sufficient 
for use in evaluating interactions with other interaction groups some distance away from the 
first group. This permits the processing of interaction data more quickly even while the data 
5 remains in a compressed format. The ability to perform numerical operations using 
compressed data allows fast processing of data using multilevel and recursive methods, as 
well as using single-level methods. 

Description of the Figures 
The advantages and features of the disclosed invention will readily be appreciated by 
10 persons skilled in the art from the following detailed description when read in conjunction 
with the drawings listed below. 
O Figure 1A illustrates a wire or rod having a physical property (e.g., a current, a 

ffl temperature, a vibration, stress, etc.) 1(1) along its length, where the shape of 1(1) is unknown. 
fk Figure IB illustrates the wire from Figure 1 A, broken up into four segments, where the 

2f 15 function 1(1) has been approximated by three known basis functions fi(l) 9 and where each 
Si basis function is multiplied by an unknown constant //. 

f%. Figure 1C illustrates a piecewise linear approximation to the function 1(1) after the 

constants I\ have been determined. 
y3 Figure 2 is a flowchart showing the process steps used to generate a compressed 

q 20 (block sparse) interaction matrix. 

Figure 3 illustrates partitioning a body into regions. 

Figure 4 shows an example of an interaction matrix (before transformation) for a body 
partitioned into five differently sized regions. 

Figure 5 shows an example of an interaction matrix after transformation (but before 
25 reordering) for a body partitioned into five regions of uniform size. 

Figure 6 shows an example of an interaction matrix after transformation and 
reordering for a body partitioned into five regions of uniform size. 

Figure 7 illustrates the block diagonal matrix D R . 
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Figure 8 is a plot showing the digits of accuracy obtained after truncating the basis 
functions for a block of the entire interaction matrix, with a block size of 67 by 93. 

Figure 9 is a plot showing the digits of accuracy obtained after truncating the basis 
functions for a block of the entire interaction matrix, with a block size of 483 by 487. 

Figure 10, consisting of Figures 10A and 10B, is a flowchart showing the process of 
generating a compressed (block sparse) impedance matrix in connection with a conventional 
moment-method computer program. 

Figure 1 1 is a three-dimensional plot showing magnitudes of the entries in a 67 by 93 
element block of the interaction matrix (before transformation) for a wire grid model using the 
method of moments. 

Figure 12 is a three-dimensional plot showing magnitudes of the entries of the 
interaction matrix from Figure 1 1 after transformation. 

In the drawings, the first digit of any three-digit number generally indicates the 
number of the figure in which the element first appears. Where four-digit reference numbers 
are used, the first two digits indicate the figure number. 

Detailed Description 

Many physical phenomena involve sources that generate a disturbance, such as an 
electromagnetic field, electromagnetic wave, a sound wave, vibration, a static field (e.g., 
electrostatic field, magnetostatic field, gravity field, etc) and the like. Examples of sources 
include a moving object (such as a loudspeaker that excites sound waves in air) and an 
electrical current (that excites electric and magnetic fields), etc. For example, the electric 
currents moving on an antenna produce electromagnetic waves. Many sources produce 
disturbances both near the source and at a distance from the source. 

Sometimes it is convenient to consider disturbances as being created by an equivalent 
source (e.g., a fictitious source) rather than a real physical source. For example, in most 
regions of space (a volume of matter for example) there are a large number of positive electric 
charges and a large number of negative electric charges. These positive and negative charges 
nearly exactly cancel each other out. It is customary to perform calculations using a fictitious 
charge, which is the net difference between the positive and negative charge, averaged over 



the region of space. This fictitious charge usually cannot be identified with any specific 
positive or negative particle. 

A magnetic current is another example of a fictitious source that is often used. It is 
generally assumed that magnetic monopoles and magnetic currents do not exist (while electric 
5 monopoles and electric currents do exist). Nevertheless, it is known how to mathematically 
relate electric currents to equivalent magnetic currents to produce the same electromagnetic 
waves. The use of magnetic sources is widely accepted, and has proven very useful for certain 
types of calculations. Sometimes, it is convenient to use a source that is a particular 
combination of electric and magnetic sources. A distribution of sources over some region of 
10 space can also be used as a source. The terms "sources" and "physical sources" are used 
herein to include all types of actual and/or fictitious sources. 
O A physical source at one location typically produces a disturbance that propagates to a 

§1 sensor (or tester) at another location. Mathematically, the interaction between a source and a 
^ tester is often expressed as a coupling coefficient (usually as a complex number having a real 
J?! 15 part and an imaginary part). The coupling coefficients between a number of sources and a 
S! number of testers is usually expressed as an array (or matrix) of complex numbers, 
fv Embodiments of this invention includes efficient methods for the computation of these 
J! complex numbers, for the storing of these complex numbers, and for computations using these 
d3 complex numbers. 

q 20 The so-called Method of Moments (MoM) is an example of numerical analysis 

procedure that uses interactions between source functions and testing functions to numerically 
solve a problem that involves finding an unknown function (that is, where the solution 
requires the determination of a function of one or more variables). The MoM is used herein 
by way of example and not as a limitation. One skilled in the art will recognize that the MoM 
25 is one of many types of numerical techniques used to solve problems, such as differential 
equations and integral equations, where one of the unknowns is a function. The MoM is an 
example of a class of solution techniques wherein a more difficult or unsolvable problem is 
broken up into one or more interrelated but simpler problems. Another example of this class 
of solution techniques is Nystrom's method. The simpler problems are solved, in view of the 
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known interrelations between the simpler problems, and the solutions are combined to 
produce an approximate solution to the original, more difficult, problem. 

For example, Figure 1A shows a wire or rod 100 having a physical property (e.g., a 
current, a temperature, a stress, a voltage, a vibration, a displacement, etc.) along its length. 
An expression for the physical property is shown as an unknown function 1(1). The problem 
is to calculate 1(1) using the MoM or a similar "divide and conquer" type of technique. By 
way of example, in many physical problems involving temperature, vibration, or electrical 
properties, etc. 1(1) will be described by an integral equation of the form: 

E(R)= jl(£)G(£,R)d£ 

Where G(£,R) is known everywhere and E(R) is known for certain values of R . In 
many circumstances, G(£,R) is a Green's function, based on the underlying physics of the 

problem, and the value of E(R) is known only at boundaries (because of known boundary 
conditions). The above equation is usually not easily solved because I(Z) is not known, and 
thus the integration cannot be performed. The above integral equation can be turned into a 
differential equation (by taking the derivative of both sides), but that will not directly provide 
a solution. Regardless of whether the above equation is expressed as an integral equation or a 
differential equation, the equation can be numerically solved for 1(1) by creating a set of 
simpler but interrelated problems as described below (provided that G(£,R) possesses certain 
mathematical properties known to those of skill in the art). 

As shown in Figure IB, in order to compute a numerical approximation for I(l\ the 
wire 100 is first divided up into four segments 101-104, and basis function f l (iXf 2 (Z), and 
f 3 (Z) are selected. In Figure IB the basis functions are shown as triangular-shaped functions 
that extend over pairs of segments. The unknown function 1(1) can then be approximated as: 

where I l9 I 2 , and 7 3 are unknown complex constants. Approximating 1(1) in this 
manner transforms the original problem from one of finding an unknown function, to a 
problem of finding three unknown constants. The above approximation for 1(1) is inserted 
into the original integral equation above to yield: 



E(R)= jl l f l (l)G(£,R)d£+ jl 2 f 2 (£)G(£,R)d£ + jl 3 f 3 (£)G(£,R)d£ 
The above integrals can now be performed because the functional form of the 
integrands are all known (G(£ 9 R) was determined by the problem being solved, the functions 
f { {l) were selected, and the constants I v I 2 and I 3 can be moved outside the integrals). 
5 However, this does not yet solve the problem because the values of I l9 I 2 and I 3 are still 
unknown. 

Fortunately, as indicated above, the value of E(R) is usually known at various 
specific locations (e.g., at boundaries). Thus, three equations can be written by selecting three 
locations i? 2 , R 3 , where the value of E(R) is known. Using these three selected 
10 locations, the above equation can be written three times as follows: 
| E(R X ) = I l \f l (£)G(£,R l )d£ + I 2 \f 2 (£)G{£j x )d£ + I 3 \f 3 {£)G(£j x )d£ 

i E(R 2 ) = \f x {£)G(£,R 2 )d£ + I 2 jf 2 (£)G(£, R 2 )d£ + I 3 jf 3 (£)G(£, R 2 )d£ 

m E(R 3 ) = I l \f 1 (£)G(£ 9 R 3 )d£ + I 2 jf 2 (£)G(£,R 3 )d£ + 1 3 jf 3 (£)G(£,R 3 )d£ 

^ Rather than selecting three specific locations for E(R) , it is known that the accuracy 

^ of the solution is often improved by integrating known values of E(R ) using a weighting 
w function over the region of integration. For example, assuming that E(R) is known along the 
y-J 15 surface of the wire 100, then choosing three weighting functions g x {i\ g 2 {^\ and g 2 (Q? the 

desired three equations in three unknowns can be written as follows (by multiplying both 

sides of the equation by gi(t) and integrating): 

\E{V) gl (e)de= i x ljMe)g x (£')G(ej')d£d£'+i 2 W^g^nc^ndide 

+i i l\f 3 (£)g 1 (nG(£,nd£de 

\E{£* )g 2 (jP )dF = /J jf x (£)g 2 {£> )G(£, £' )d£d£'+I 2 { J/ 2 (£)g 2 (£ f )G(£ 9 t )dftW 

+ I 3 \\f 3 (£)g 2 {£)G{£,nd£d£ 
lEing^W^I, \ jM£)g 3 (P)G(£J')d£d£<+I 2 jlf 2 (£)g 3 (nG(£,f)d£df 

+i 3 \\M£)g 3 (nG(£,nd£dt 
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Note that the above double-integral equations reduce to the single-integral forms if the 
weighting functions gj(l) are replaced with delta functions. 

The three equations in three unknowns can be expressed in matrix form as: 

V=ZI 

or 















v 2 










h 


r 3 _ 




. Z 3. 


Z 32 


^33 _ 





where 
and 

Solving the matrix equation yields the values of I l9 7 2? and 7 3 . The values /„ 7 2 , and I 3 
can then be inserted into the equation I(£)& A/W + ^W + ^/sW to g* ve 311 
approximation for /(£). If the basis functions are triangular functions as shown in Figure 1B ? 
then the resulting approximation for 1(1) is a piecewise linear approximation as shown in 
Figure 1C. The If are the unknowns and the Vf are the conditions (tpyically, the Vj are 
knowns). Often there are the same number of conditions as unknowns. In other cases, there 
are more conditions than unknowns or less conditions than unknown. 

The accuracy of the solution is largely determined by the shape of the basis functions, 
by the shape of the weighting functions, and by the number of unknowns (the number of 
unknowns usually corresponds to the number of basis functions). 

Unlike the Moment Method described above, some techniques do not use explicit 
basis functions, but, rather, use implicit basis functions or basis-like functions. For example, 
Nystrom's method produces a numerical value for an integral using values of the integrand at 
discrete points and a quadrature rule. Although Nystrom's method does not explicitly use an 
expansion in terms of explicit basis functions, nevertheless, in a physical sense, basis 
functions are still being used (even if the use is implicit). That is, the excitation of one 



unknown produces some reaction throughout space. Even if the computational method does 
not explicitly use a basis function, there is some physical excitation that produces 
approximately the same reactions. All of these techniques are similar, and one skilled in the 
art will recognize that such techniques can be used with the present invention. Accordingly, 
the term "basis function" will be used herein to include such implicitly used basis functions. 
Similarly, the testers may be implicitly used. 

When solving most physical problems (e.g., current, voltage, temperature, vibration, 
force, etc), the basis functions tend to be mathematical descriptions of the source of some 
physical disturbance. Thus, the term "source" is often used to refer to a basis function. 
Similarly, in physical problems, the weighting functions are often associated with a receiver 
or sensor of the disturbance, and, thus, the term "tester" is often used to refer to the weighting 
functions. 

As described above in connection with Figures 1A-1C, in numerical solutions, it is 
often convenient to partition a physical structure or a volume of space into a number of 
smaller pieces and associate the pieces with one or more sources and testers. In one 
embodiment, it is also convenient to partition the structure of (or volume) into regions, where 
each region contains a group of the smaller pieces. Within a given region, some number of 
sources is chosen to describe with sufficient detail local interactions between sources and 
testers within that region. A similar or somewhat smaller number of sources in a given region 
is generally sufficient to describe interactions between sources in the source region and testers 
in the regions relatively close by. When the appropriate sources are used, an even smaller 
number of sources is often sufficient to describe interactions between the source region and 
testers in regions that are not relatively close by (i.e., regions that are relatively far from the 
source region). 

Embodiments of the present invention include methods and techniques for finding 
composite sources. Composite sources are used in place of the original sources in a region 
such that a reduced number of composite sources is needed to calculate the interactions with a 
desired accuracy. 
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In one embodiment, the composite sources for a first region are the same regardless of 
whether the composite sources in the first region are interacting with a second region, a third 
region, or other regions. The use of the same composite sources throughout leads to efficient 
methods for factoring and solving the interaction matrix. 

Considering the sources in the first region, one type of source is the so-called 
multipole, as used in a multipole expansion. Sources like wavelets are also useful. In some 
cases wavelets allow a reduced number of composite sources to be used to describe 
interactions with distant regions. However, there are disadvantages to wavelet and multipole 
approaches. Wavelets are often difficult to use, and their use often requires extensive 
modifications to existing or proposed computer programs. Wavelets are difficult to implement 
on non-smooth and non-planar bodies. 

Multipole expansions have stability problems for slender regions. Also, while a 
multipole expansion can be used for describing interactions with remote regions, there are 
severe problems with using multipoles for describing interactions within a region or between 
spatially close regions. This makes a factorization of the interaction matrix difficult. It can be 
very difficult to determine how to translate information in an interaction matrix into a wavelet 
or multipole representation. 

Figure 2 is a flowchart that illustrates a compression technique 200 for compressing an 
interaction matrix by combining groups of sources and groups of testers into composite 
sources and testers. The use of composite sources and composite testers allows the original 
interaction matrix to be transformed into a block sparse matrix having certain desirable 
properties. 

Embodiments of the present invention include a technique for computing and using 
composite sources to provide compression of an interaction matrix by transforming the 
interaction matrix into a block sparse matrix. The present technique is compatible with 
existing and proposed computer programs. It works well even for rough surfaces and irregular 
grids of locations. For a given region, the composite sources allow computation of a 
disturbance (e.g., radiation) produced by the source throughout a desired volume of space. A 
reduced number of these composite sources is sufficient to calculate (with a desired accuracy) 



disturbances at other relatively distant regions. This method of compressing interaction data 
can be used with a variety of computational methods, such as, for example, an LU (Lower 
Triangular Upper triangular) factorization of a matrix or as a preconditioned conjugate 
gradient iteration. In many cases, the computations can be done while using the compressed 
storage format. 

Figure 2 is a flowchart 200 illustrating the steps of solving a numerical problem using 
composite sources. The flowchart 200 begins in a step 201 where a number of original 
sources and original testers are collected into groups, each group corresponding to a region. 
Each element of the interaction matrix describes an interaction (a coupling) between a source 
and a tester. The source and tester are usually defined, in part, by their locations in space. 
The sources and testers are grouped according to their locations in space. In one embodiment, 
a number of regions of space are defined. A reference point is chosen for each region. 
Typically the reference point will lie near the center of the region. The sources and testers are 
grouped into the regions by comparing the location of the source or tester to the reference 
point for each region. Each source or tester is considered to be in the region associated with 
the reference point closest to the location. (For convenience, the term "location" is used 
hereinafter to refer to the location of a source or a tester.) 

Other methods for grouping the sources and testers (that is, associating locations with 
regions) can also be used. The process of defining the regions is problem-dependent, and in 
some cases the problem itself will suggest a suitable set of regions. For example, if the 
sources and testers are located on the surface of a sphere, then curvilinear-square regions are 
suggested. If the sources and testers are located in a volume of space, then cubic regions are 
often useful. If the sources and testers are located on a complex three-dimensional surface, 
then triangular patch-type regions are often useful. 

Generally the way in which the regions are defined is not critical, and the process used 
to define the regions will be based largely on convenience. However, it is usually preferable 
to define the regions such that the locations of any region are relatively close to each other, 
and such that there are relatively few locations from other regions close to a given region. In 
other words, efficiency of the compression algorithm is generally improved if the regions are 
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as isolated from one another as reasonably possible. Of course, adjacent regions are often 
unavoidable, and when regions are adjacent to one another, locations near the edge of one 
region will also be close to some locations in an adjacent region. Nevertheless, the 
compression will generally be improved if, to the extent reasonably possible, regions are 
defined such that they are not slender, intertwining, or adjacent to one another. For example, 
Figure 3 illustrates a volume of space partitioned into a rectangular box 300 having eleven 
regions A through K corresponding to reference points 301-31 1. 

As shown in Figure 2, after the step 201 the process advances to a step 202. In the 
step 202, the unknowns are renumbered, either explicitly or implicitly, so that locations within 
the same region are numbered consecutively. It is simpler to continue this description as if the 
renumbering has actually been done explicitly. However, the following analysis can also be 
performed without explicit renumbering. 

The term "spherical angles" is used herein to denote these angles. One skilled in the 
art will recognize that if a two-dimensional problem is being solved, then the spherical angles 
reduces to a planar angle. Similarly, one skilled in the art will recognize that if a higher- 
dimensional problem is being solved (such as, for example, a four dimensional space having 
three dimensions for position and one dimension for time) then the term spherical angle 
denotes the generalization of the three-dimensional angle into four-dimensional space. Thus, 
in general, the term spherical angle is used herein to denote the notion of a "space-filling" 
angle for the physical problem being solved. 

After renumbering, the process advances to a block 203 where one or more composite 
sources for each region are determined. If there are p independent sources within a region, 
then q composite sources can be constructed (where q<p). The construction of composite 
sources begins by determining a relatively dense set of far-field patterns (usually described in 
a spherical coordinate system) at relatively large distances from the region. As used herein, 
far-field refers to the field in a region where the field can be approximated in terms of an 
asymptotic behavior. For example, in one embodiment, the far-field of an antenna or other 
electromagnetic radiator includes the field at some distance from the antenna, where the 
distance is relatively larger than the electrical size of the antenna. 
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A far-field pattern is constructed for each independent source. In the present context, 
dense means to avoid having any overly-large gaps in the spherical angles used to calculate 
the set of disturbances. Dense also means that if the disturbance is represented by a vector, 
then each vector component is represented. For example, for a scalar problem, one might 
choose p spherical angles. These angles are typically substantially equally spaced, and the 
ranges of angles include the interaction angles occurring in the original interaction matrix (if 
all of the interactions described in the original matrix lie within a plane, then one can choose 
directions only within that plane rather than over a complete sphere). 

The far-field data is stored in a matrix s having p columns (one column for each source 
location within the region), and rows associated with angles. While each source is logically 
associated with a location in a given region, these sources are not necessarily located entirely 
within that region. While each source corresponds to a location (and each location is assigned 
to a region), sources that have a physical extent can extend over more than one region. The 
entries in the matrix s can be, for example, the field quantity or quantities that emanate from 
each source. It is desirable that the field quantity is chosen such that when it (or they) are zero 
at some angle then, to a desired approximation, all radiated quantities are zero at that angle. 
While it is typically desirable that the angles be relatively equally spaced, large deviations 
from equal spacing can be acceptable. 

These composite sources are in the nature of equivalent sources. A smaller number of 
composite sources, compared to the number of sources they replace, can produce similar 
disturbances for regions of space removed from the region occupied by these sources. 

As described above, sources are collected into groups of sources, each group being 
associated with a region. For each group of sources, a group of composite sources is 
calculated. The composite source is in the nature of an equivalent source that, in regions of 
space removed from the region occupied by the group in replaces, produces a far-field 
(disturbance) similar to the field produced by the group it replaces. Thus, a composite source 
(or combination of composite sources) efficiently produces the same approximate effects as 
the group of original sources at desired spherical angles and at a relatively large distance. To 
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achieve a relatively large distance, is it often useful to use a limiting form as the disturbance 
goes relatively far from its source. 

Each composite source is typically a linear combination of one or more of the original 
sources. A matrix method is used to find composite sources that broadcast strongly and to 
5 find composite sources that broadcast weakly. These composite sources are constructed from 
the original sources. The matrix method used to find composite sources can be a rank- 
revealing factorization such as singular value decomposition. For a singular value 
decomposition, the unitary transformation associated with the sources gives the composite 
sources as a linear combination of sources. 

10 Variations of the above are possible. For example, one can apply the singular value 

decomposition to the transpose of the s matrix. One can employ a Lanczos Bi-diagonalization, 
or related matrix methods, rather than a singular value decomposition. There are other known 
methods for computing a low rank approximation to a matrix. Some examples of the use of 
Lanczos Bidiagonalization are given in Francis Canning and Kevin Rogovin, "Fast Direct 

15 Solution of Standard Moment-Method Matrices," IEEE AP Magazine, Vol. 40, No. 3, June 
1998, pp. 15-26. 

There are many known methods for computing a reduced rank approximation to a 
matrix. A reduced rank approximation to a matrix is also a matrix. A reduced rank matrix 
with m columns can be multiplied by any vector of length m. Composite sources that 
20 broadcast weakly are generally associated with the space of vectors for which that product is 
relatively small (e.g., in one embodiment, the product is zero or close to zero). Composite 
sources that broadcast strongly are generally associated with the space of vectors for which 
that product is not necessarily small. 

Composite sources can extend over more than one region. In one embodiment, this is 
25 achieved by using the technique used with Malvar wavelets (also called local cosines) to 
extend Fourier transforms on disjoint intervals to overlapping orthogonal functions. 

Persons of ordinary skill in the art know how near-field results are related to far-field 
results. A relationship between near-field and far-field can be used in a straightforward way to 
transform the method described above using far-field data into a method using near- field data. 
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Note that, the "far-field" as used herein is not required to correspond to the traditional 2cf/X 
far-field approximation. Distances closer than 2cf/X can be used (although closer distances 
will typically need more composite sources to achieve a desired accuracy). A distance 
corresponding to the distance to other physical regions is usually far enough, and even shorter 
distances can be acceptable. 

Once composite sources are found, the process advances to a step 204 where 
composite testers are found. Composite testers are found in a manner analogous to the way 
that composite sources are found. Recall that composite sources are found using the way in 
which sources of the interaction matrix "broadcast" to distant locations. Composite testers are 
found using the way in which the testers of the interaction matrix "receive" from a dense 
group of directions for a distant disturbance. It is helpful if the received quantity or quantities 
which are used include relatively all field quantities, except (optionally) those which are very 
weakly received. For example, when receiving electromagnetic radiation from a distant 
source, the longitudinal component is approximately zero and can often be neglected. A 
matrix R describing how these testers receive is formed. A matrix method is used to construct 
composite testers that receive strongly and testers that receive weakly. The matrix method 
can be a rank-revealing factorization such as singular value decomposition. A singular value 
decomposition gives the composite testers as a linear combination of the testers which had 
been used in the original matrix description. 

Once composite sources and testers have been found, the process advances to a step 
205 or to a step 215 where the interaction matrix is transformed to use composite sources and 
testers. The steps 205 and 215 are alternatives. Figure 4 shows an example of an interaction 
matrix 400 having 28 unknowns (28 sources and 28 testers) grouped into five physical regions 
(labeled I-V). The shaded block 401 of the matrix 400 represents the interaction for sources in 
the fourth region (region IV) and testers in the second region (region II). The interaction of a 
pair of regions describes a block in the interaction matrix 400. The blocks of the transformed 
matrix can be computed at any time after the composite functions for their source and tester 
regions are both found. That is, the block 401 can be computed after composite sources for 
region IV and testers for region II are found. 
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The step 215 of Figure 2 shows one method for computing all of the blocks in the 
matrix 400 by computing the entries for these blocks using the original sources and testers. 
Then, the process advances to an optional step 216 where these blocks are transformed into a 
description in terms of the composite sources and composite testers. 

One advantage of using composite sources and testers is that many entries in the 
transformed matrix will be zero. Therefore, rather than transforming into a description using 
composite modes, the step 205 shows calculation of the transformed block directly using the 
composite sources and composite testers (without first calculating the block using the original 
sources and testers). In other words, the composite sources are used as basis functions, and 
the composite testers are used as weighting functions. Within each block, entries that are 
known au priori to be zero (or very small) are not calculated. 

Further savings in the storage required are possible. After each block has been 
transformed, only the largest elements are kept. No storage needs to be used for the elements 
that are approximately zero. Many types of block structures, including irregular blocks and 
multilevel structures, can also be improved by the use of this method for storing a block 
sparse matrix. This will usually result in a less regular block structure. As an alternative, it is 
also possible to store a portion of the interaction data using composite sources and testers and 
to store one or more other portions of the data using another method. 

The non-zero elements of the interaction matrix typically occur in patterns. After 
either the step 205 or the step 216, the process advances to a step 206 where the interaction 
matrix is reordered to form regular patterns. For a more uniform case, where all of the regions 
have the same number of sources, the resulting transformed matrix T is shown in Figure 5. 
Figure 5 shows non-zero elements as shaded and zero elements as unshaded. If only a 
compressed storage scheme is desired, the process can stop here. However, if it is desired to 
calculate the inverse of this matrix, or something like its LU (lower-upper triangular) 
factorization, then a reordering can be useful. 

The rows and columns of the interaction matrix can be reordered, to produce a matrix 
T A in the form shown in Figure 6. This permutation moves the composite sources that 
broadcast strongly to the bottom of the matrix, and it moves the composite testers which 
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receive strongly to the right side of the matrix. The interaction between composite sources 
and composite testers is such that the sizes of the matrix elements can be estimated au priori. 
A matrix element that corresponds to an interaction between a composite source that radiates 
strongly and a composite tester that receives strongly will be relatively large. A matrix 
element that corresponds to an interaction between a composite source that radiates strongly 
and a composite tester that receives weakly will be relatively small. Similarly, a matrix 
element that corresponds to an interaction between a composite source that radiates weakly 
and a composite tester that receives strongly will be relatively small. A matrix element that 
corresponds to an interaction between a composite source that radiates weakly and a 
composite tester that receives weakly will be very small. 

The permuted matrix T A often will tend to be of a banded form. That is, the non-zero 
elements down most of the matrix will tend to be in a band near the diagonal. For a matrix of 
this form, there are many existing sparse-matrix LU factorers and other matrix solvers, that 
work well. The order shown in Figure 6 is one example. The LU decomposition of the matrix 
T A can be computed very rapidly by standard sparse matrix solvers. The matrices L and U of 
the LU decomposition will themselves be sparse. For problems involving certain types of 
excitations, only a part of the matrices L and U will be needed, and this can result in further 
savings in the storage required. 

After reordering, the process 200 advances to a step 207 where the linear matrix 
problem is solved. The matrix problem to be solved is written as: 

T A G=H 

where the vector H represents the excitation and the vector G is the desired solution for 
composite sources. The excitation is the physical cause of the sound, temperature, 
electromagnetic waves, or whatever phenomenon is being computed. If the excitation is very 
distant (for example, as for a plane wave source), H will have a special form. If the vector H 
is placed vertically (as a column vector) alongside the matrix of Figure 6, the bottom few 
elements of H alongside block 602, will be relatively large, and the remaining elements of H 
will be approximately equal to zero. The remaining elements of H are approximately zero 
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because the composite testers separate the degrees of freedom according to how strongly they 
interact with a distant source. 

When T A is factored by LU decomposition, then: 

T A = LU; 
LUG = H; 

and this is solved by the following two-step process; 

Step I: Find X in LX = H 

Step II: Find Gin UG = X 

The matrix L is a lower triangular matrix (meaning elements above its diagonal are 
zero). It follows immediately from this that if only the bottom few elements of H are non- 
zero, then only the bottom elements of X are non-zero. As a consequence, only the bottom 
right portion of L is needed to compute G. The remaining parts of L were used in computing 
this bottom right portion, but need not be kept throughout the entire process of computing the 
LU decomposition. This not only results in reduced storage, but also results in a faster 
computation for Step I above. 

If only the far-field scattered by an object needs to be found, then further efficiencies 
are possible. In that case, it is only necessary to find the bottom elements of G, corresponding 
to the bottom non-zero elements of H. This can be done using only the bottom right portion 
of the upper triangular matrix U. This results in efficiencies similar to those obtained for L. 

For other types of excitations, similar savings are also possible. For example, for 
many types of antennas, whether acoustic or electromagnetic, the excitation is localized 
within one active region, and the rest of the antenna acts as a passive scattered In that case, 
the active region can be arranged to be represented in the matrix of Figure 6 as far down and 
as far to the right as possible. This provides efficiencies similar to those for the distant 
excitation. 

A permutation of rows and a permutation of columns of the matrix T of Figure 5 
brings it to the matrix T A of Figure 6. These permutations are equivalent to an additional 
reordering of the unknowns. Thus, a solution or LU decomposition of the matrix T A of Figure 
6 does not immediately provide a solution to the problem for the original data. Several 
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additional steps are used. These steps involve: including the effects of two permutations of the 
unknowns; and also including the effect of the transformation from the original sources and 
testers to using the composite sources and composite testers. 

Direct methods (such as LU decomposition) and iterative methods can both be used to 
solve the matrix equation herein. An iterative solution, with the compressed form of the 
matrix, can also be used with fewer computer operations than in the prior art. Many iterative 
methods require the calculation of the product of a matrix and a vector for each iteration. 
Since the compressed matrix has many zero elements (or elements which may be 
approximated by zero), this can be done more quickly using the compressed matrix. Thus, 
each iteration can be performed more quickly, and with less storage, than if the uncompressed 
matrix were used. 

The compressed format of T A has an additional advantage. In many cases, there is a 
way to substantially reduce the number of iterations required, resulting in further increases in 
speed. For example, in the method of conjugate gradients, the number of iterations required 
to achieve a given accuracy depends on the condition number of the matrix. (The condition 
number of a matrix is defined as its largest singular value divided by its smallest.) Physical 
problems have a length scale, and one interpretation of these composite sources and 
composite testers involves length scales. These composite sources and composite testers can 
be described in terms of a length scale based on a Fourier transform. This physical fact can be 
used to improve the condition number of the matrix and therefore also improve the speed of 
convergence of the iterative method. 

A composite source is a function of spatial position, and its Fourier transform is a 
function of "spatial frequency." Composite sources that broadcast weakly tend to have a 
Fourier transform that is large when the absolute value of this spatial frequency is large. 
There is a correlation between how large this spatial frequency is and the smallness of the 
small singular values of the matrix. This correlation is used in the present invention to 
provide a method to achieve convergence in fewer iterations. 

Two matrices, P R and P L are defined as right and left preconditioning matrices to the 
compressed matrix. Each column of the compressed matrix is associated with a composite 



-20- 



source. This composite source can be found using a matrix algebra method, such as a rank- 
revealing factorization (e.g., singular value decomposition and the like). The rank-revealing 
factorization method provides some indication of the strength of the interaction between that 
composite source and other disturbances. For example, using a singular value decomposition, 
the associated singular value is proportional to this strength. The diagonal matrix P R is 
constructed by choosing each diagonal element according to the interaction strength for the 
corresponding composite source. The diagonal element can be chosen to be the inverse of the 
square root of that strength. Similarly, the diagonal matrix P L can be constructed by choosing 
each diagonal element according to the interaction strength for its associated composite tester. 
For example, the diagonal element can be chosen to be the inverse of the square root of that 
strength. 

If the compressed matrix is called T, then the preconditioned matrix is 

p = p L 'j 1 p R 

The matrix P will often have a better (i.e., smaller) condition number than the matrix 
T. There are many iterative methods that will converge more rapidly when applied to the 
preconditioned matrix P rather than to T. 

One embodiment of the composite source compression technique is used in connection 
with the computer program NEC2. This program was written at Lawrence Livermore National 
Laboratory during the 1970s and early 1980s. The NEC2 computer program itself and 
manuals describing its theory and use are freely available over the Internet. The following 
development assumes NEC2 is being used to calculate the electromagnetic fields on a body 
constructed as a wire grid. 

NEC2 uses electric currents flowing on a grid of wires to model electromagnetic 
scattering and antenna problems. In its standard use, NEC2 generates an interaction matrix, 
herein called the Z matrix. The actual sources used are somewhat complicated. There is at 
least one source associated with each wire segment. However, there is overlap so that one 
source represents current flowing on more than one wire segment. NEC2 uses an array CURX 
to store values of the excitation of each source. 
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Figure 10 is a flowchart 1000 showing the process of using NEC2 with composite 
sources and composite testers. The flowchart 1000 begins at a step 1001 where the NEC2 
user begins, as usual, by setting up information on the grid of wires and wire segments. The 
process then advances to a step 1002 to obtain from NEC2 the number of wire segments, their 
5 locations (x,y,z coordinates), and a unit vector i for each segment. The vector t is tangent 
along the wire segment, in the direction of the electric current flow on the wire segment. 
' Next, in a step 1003, the wire grid is partitioned into numbered regions. A number of 

reference points are chosen. The reference points are roughly equally spaced over the volume 
occupied by the wire grid. Each wire segment is closest to one of these reference points, and 
10 the segment is considered to be in the region defined by the closest reference point. In one 
embodiment, the number of such points (and associated regions) is chosen as the integer 
y closest to the square root of N (where N is the total number of segments). This is often an 
CP effective choice, although the optimum number of points (and associated regions) depends on 
m many factors, and thus other values can also be used. For a set of N segments, each wire 
~1 15 segment has an index, running from 1 to N. 

H After the step 1003, the process advances to a step 1004 where the wires are sorted by 

P region number. After sorting, the numbering of the wires is different from the numbering used 
by NEC2. Mapping between the two numbering systems is stored in a permutation table that 
translates between these different indexes for the wire segments. Using this new numbering 
□ 20 of segments, an array "a" is created, such that a(p) is the index of the last wire segment of the 
p th region (define a(0) = 0 for convenience). 

After renumbering, the process advances to a step 1005 where composite sources and 
composite testers for all regions are calculated. Source region p corresponds to unknowns 
a(p-l) + 1 through a(p) in the ordering. Define MasM= a(p) - a(p - 1). Choose M directions 
25 substantially equally spaced throughout three-dimensional space. In other words, place M 
roughly equally spaced points on a sphere, and then consider the M directions from the center 
of the sphere to each point. The order of the directions is unimportant. One convenient method 
for choosing these points is similar to choosing points on the earth. For example, choose the 
North and South poles as points. A number of latitudes are used for the rest of the points. For 
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each chosen latitude, choose points equally spaced at a number of longitudes. This is done so 
that the distance along the earth between points along a latitude is approximately the same as 
the distance between the latitude lines holding the points. It is desirable that the points are 
equally spaced. However, even fairly large deviations from equal spacing are tolerable. 

Now generate a matrix A of complex numbers with 2M rows and M columns. For 
m = 1 to M and for n =1 to M, compute elements of this matrix two at a time: the element at 
row m and column n and also the element at row m + M and column n. To compute these two 
elements, first fill the NEC2 array CURX with zero in every position. Then, set position a(p- 
1) + n of CURX to unity. A value of unity indicates that only source number a(p-l) + n is 
excited. This source is associated with the wire segment of that number, even though it 
extends onto neighboring segments. The matrix Z is defined in terms of these same sources. 
Then, call the NEC2 subroutine CABC (CURX) . The subroutine CABC generates a different 
representation of the source, but the same representation that the NEC2 subroutine FFLD 
uses. This representation is automatically stored within NEC2. The m th direction previously 
chosen can be described in spherical coordinates by the pair of numbers (Theta, Phi) . 
After calculating Theta and Phi, the NEC2 subroutine FFLD (Theta, Phi, ETH, EPH) is 
called. Theta and Phi are inputs, as are the results from CABC. The outputs from FFLD are 
the complex numbers ETH and EPH. ETH and EPH are proportional to the strengths of the 
electric field in the far-field (far away from the source) in the theta and phi directions 
respectively. ETH is placed in row m and column n 9 (m,ri) 9 of A. EPH is placed at row m+M 
and column n of A. Alternatively, there are other ways to compute the numbers ETH and EPH 
produced by FFLD. For example, it will apparent to one of ordinary skill in the art that the 
subroutine FFLD can be modified to produce an answer more quickly by making use of the 
special form of the current, since most of the entries in the current are zero. 

Next, a singular value decomposition of A is performed, such that: 

A = UDV h 

where U and V are unitary matrices, and D is a diagonal matrix. The matrix U will not 
be used, so one can save on computer operations by not actually calculating IL The matrix V 
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has M rows and M columns. Since these calculations were performed for the p region, the 
square matrix dp R is defined by 

d*= V 

The reason for this choice comes from the fact that 

AV = UD 

and that each successive columns of the product UD tends to become smaller in magnitude. 
They become smaller because U is unitary and the singular values on the diagonal of D 
decrease going down the diagonal. 

Next, assemble an N by N block diagonal matrix D R . That is, along the diagonal the 
first block corresponds to p=l. Starting at the bottom right corner of that block, attach the 
block for p-2 9 etc., as shown in Figure 7. 

Next a similar procedure is followed to find the block diagonal matrix D L . For each 
region p, a matrix A is filled as before. However, this time this region is considered as 
receiving rather than as transmitting. Again A will have 2Mrows and M columns, where M = 
a(p) - a(p-l). Again there are M directions, but now those are considered to be the receiving 
directions. 

To understand what is to be put into A, it is instructive to note how the NEC2 
computer program defines the interaction matrix Z. When used with wire grid models, the 
sources radiate electric and magnetic fields. However, it is the electric field reaching another 
segment that is used in NEC2. Each matrix element of Z is computed by computing the 
component of that electric field which is in the direction of the tangent to the wire segment. 

For the pair of numbers (m,n) 9 where m = 1,..., M and n = 1,...,M, the matrix entries for 

A at (m,n) and {m+M,ri) are calculated as follows. Compute a unit vector k in the 

direction. Find the unit vector tangent to segment number n, and call it i . The position of the 
center of wire segment number n is found and is designated as the vector X. Then compute 
the vector 

f = (t-(k't)k)e J2 ^ u 
where the wavelength is given by X (NEC2 uses units where X = 1). 
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Note that the Green's function for this problem has a minus sign in the exponential, 

and the foregoing expression does not. This is because the direction of k is outward, which is 
opposite to the direction of propagation of the radiation. 

For problems in electromagnetics, the physical wavelength X is greater than zero. If a 
problem in electrostatics were being solved instead, electrostatics can be considered as the 
limit when the wavelength becomes arbitrarily large. The complex exponential above can 
then be replaced by unity. Also, for electrostatics, the relevant field quantity can be 

longitudinal (meaning f would be parallel to k ). 

For this value of m (and associated direction k% spherical coordinates define two 
directions called the theta and the phi directions. These directions are both perpendicular to 

the direction of k . Compute the components of f in each of these directions, and designate 
them as f Theta and f Phi. These are complex numbers. Then place f Theta in row m and 
column n of A and place f Phi in row m + Mand column n of A. 

The matrix A is a matrix of complex numbers. Take the complex conjugate of A, 
(A*), and perform a singular value decomposition on it, such that: 

A* = UDV h 

Now define the left diagonal block for region p, dp L , as 

dp L = V h 

The superscript h on V, indicates Hermitian conjugate. The definition of the blocks for 
the right side did not have this Hermitian conjugate. From these diagonal blocks, assemble an 
N by N matrix D L in the same way as D R was assembled earlier. The motivation for these 
choices is partly that the matrix DU h has rows that tend to become smaller. Further, it is 
expected that the Green's function that was used in creating Z has properties similar to the 
far-field form used in creating A*. The formula 

V h A 4 = DU h 

shows that V h A* will also have successive rows that tend to become smaller. The choices 
described above suggest that successive rows of each block of the compressed matrix will also 
have that property. 
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It should be noted that the matrix A, whether used for the right side or for the left side, 
can be filled in other ways as well. For example, with an appropriate (consecutive in space) 
ordering of the angles, A can be made as an M by M matrix by using theta polarization 
(f Theta) values for one angle and phi polarization values (f Phi) for the next. Usually, it is 
desirable that A does not leave large gaps in angle for any component of the electric field, 
which is important far from the source or receiver. 

In performing the singular value decompositions for the right and left sides, singular 
values are found each time. Figures 8 and 9 show the singular values found for blocks of size 
67 by 93 and 483 by 487, respectively. These calculations were done for a wire grid model 
with NEC2. The singular values are plotted in terms of how many orders of magnitude they 
are smaller than the largest singular value, and this is called "Digits of Accuracy" on the plots. 
Figures 8 and 9 show the accuracy that is achieved when truncating to a smaller number of 
composite sources or composite testers for regions that are relatively far apart. For regions 
that are closer together, the desired accuracy often requires the information from more 
composite sources and composite testers to be kept. 

After computing composite sources and composite testers, the process advances to a 
step 1006 where a new matrix T, which uses the composite sources and testers associated with 
D L and D R , is computed. The matrix is T given by the equation 

T - D L Z D R 

T can be efficiently generated by using the numbering of the wire segments developed 
herein (rather than the numbering used in NEC2). The matrix Z is computed by NEC2 and 
renumbered to use the numbering described herein. Note that a block structure has been 
overlaid on Z and T. This block structure follows from the choice of regions. Figure 4 shows 
one example of a block structure. Block {p f q} of the matrix T, to be called T{p,q}, is the part 
of T for the rows in region number p and the columns in region number q. The formula for T 
given above is such that T{p f q) only depends on Z{p,q}. Thus, only one block of Z at a time 
needs to be stored. 

In the step 1006, T is assembled one block at a time. For each block of T, first obtain 
from NEC2 the corresponding block of Z. The wire segments within a block are numbered 
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consecutively herein (NEC2 numbers them differently). Thus, first renumber Z using the 
renumber mapping from step 1004, and then perform the calculation: 

1{p,q} = l£Z{p,q}&? 

Many of the numbers in T{p,q} will be relatively small. An appropriate rule based on 
a desired accuracy is used to choose which ones may be approximated by zero. The remaining 
non-zero numbers are stored. Storage associated with the zero-valued elements of T{p,q) and 
of Z{p,q} can be released before the next block is calculated. The top left portion of T 
has matrix elements which will be kept. Anticipating this, the calculation speed can be 
increased by not calculating either the right portion or the bottom portion of T{p,q}. 

The matrix T is a sparse matrix, and it can be stored using an appropriate data 
structure for a sparse matrix. For a matrix with N z non-zero elements, an array Z z (i) for 
/ = 1,...,^ can be used, where Z z {i) is the complex value of the f h matrix element. There are 
two integer valued arrays, I z (i) and J z (i) for / = l,...,N z . I z (i) gives the row number where the 
i th matrix element occurs in T and J z (i) its column number. 

After calculation of T, the process proceeds to a process block 1007 where the rows 
and columns of the matrix T are reordered to produce a matrix T A . The matrix T is reordered 
into a matrix T A so that the top left corner of every block of T A ends up in the bottom right 
corner of the whole matrix. The T A form is more amenable to LU factorization. Figure 5 
shows an example of a matrix T, and Figure 6 shows an example of a matrix T A after 
reordering. One embodiment uses a solver that has its own reordering algorithms thus 
negating the need for an explicit reordering from T to T A . 

After reordering, the process advances to a step 1008 where the matrix T A is passed to 
a sparse matrix solver, such as, for example, the computer program "Sparse," from the 
Electrical Engineering Department of University of California at Berkeley. The program 
Sparse can be used to factor the matrix T A into a sparse LU decomposition. 

NEC2 solves the equation 

J - Z 1 E 

for various vectors E. In Figure 10, the solution of the above matrix equation is done 
in steps 1009-1016 or, alternatively, in steps 1017-1023. The sequence of steps 1009-1016 is 
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used with a matrix equation solver that does not provide reordering. The sequence of steps 
1017-1023 is used with a matrix equation solver that does provide reordering. 

In the step 1009, the vector E is computed by NEC2. Then, in the step 1010, the 
elements of E are permutated (using the same permutation as that used in the step 1004) to 
produce a vector E\ This permutation is called the region permutation. Next, in the step 
1011, E' is expressed in terms of composite testers by multiplying E' by D L , giving D L E\ 
Then, in the step 1012, the same permutation used in the step 1007 is applied to D L E' to yield 
(D L E') A . (This permutation is called the solver permutation.) Then, in the step 1013, a matrix 
equation solver (such as, for example, the solver known as "SPARSE") is used to solve for the 
vector Y A from the equation 

T A (Y A )=(D L E ' ) A 

Then, in the step 1014, the inverse of the solver permutation is applied to Y A to yield 
Y. In the subsequent step 1015, J' is computed from the equation 

J'=D R Y 

In the subsequent, and final, step 1016, the inverse of the region permutation is applied 
to J' to yield the desired answer J. 

Alternatively, the embodiment shown in steps 1017-1023 is conveniently used when 
the matrix equation solver provides its own reordering algorithms, thus eliminating the need 
to reorder from T to T A (as is done in the step 1007 above). In the step 1017, a reordering 
matrix solver is used to solve the matrix T. In the subsequent step 1018, the vector E is 
computed by NEC2. Then, in the step 1019, the elements of E are permutated using the 
region permutation to produce a vector E\ Then, in the step 1020, D L E' is computed. The 
process then proceeds to the step 1021 where the equation 

TY=D L E' 

is solved for Y. After Y is computed, the process advances to the step 1022 where J' 
is calculated from the equation 

J'=D R Y 

Finally, in the step 1023, the inverse of the region permutation is applied to J' to yield 
the desired answer J. 
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Many matrix elements are made small by this method. Figures 1 1 and 12 show before 
and after results for a problem using a wire grid model in NEC2, with a matrix Z of size 2022 
by 2022 and a block of size 67 by 93. Figure 1 1 shows the magnitudes of the matrix elements 
before changing the sources and testers, meaning it shows a 67 by 93 block of the renumbered 
Z. Figure 12 shows this same block of T. The matrix T has a regular structure wherein the 
large elements are in the top left corner. This is a general property of the transformed matrix. 
For larger blocks, the relative number of small matrix elements is even better. 

The algorithms expressed by the flowchart shown in Figure 2 can be implemented in 
software and loaded into a computer memory attached to a computer processor to calculate, 
for example, propagation of energy, pressure, vibration, electric fields, magnetic fields, strong 
nuclear forces, weak nuclear forces, etc. Similarly, the algorithms expressed by the flowchart 
shown in Figure 10 can be implemented in software and loaded into a computer memory 
attached to a computer processor to calculate, for example, electromagnetic radiation by an 
antenna, electromagnetic scattering, antenna properties, etc. 

Although the foregoing has been a description and illustration of specific 
embodiments of the invention, various modifications and changes can be made thereto by 
persons skilled in the art without departing from the scope and spirit of the invention. For 
example, in addition to electromagnetic fields, the techniques described above can also be 
used to compress interaction data for physical disturbances involving a heat flux, an electric 
field, a magnetic field, a vector potential, a pressure field, a sound wave, a particle flux, a 
weak nuclear force, a strong nuclear force, a gravity force, etc. The techniques described 
above can also be used for lattice gauge calculations, economic forecasting, state space 
reconstruction, and image processing (e.g., image formation for synthetic aperture radar, 
medical, or sonar images). Accordingly, the invention is limited only by the claims that 
follow. 
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WHAT IS CLAIMED IS : 

1 . A method of data compression, comprising: 

partitioning a first set of basis functions into groups, each group corresponding 
to a region, each basis function corresponding to one unknown in a system of linear 
equations, each of said basis functions corresponding to an original source; 

selecting a plurality of spherical angles; 

calculating a far-field disturbance produced by each of said basis functions in a 
first group for each of said spherical angles to produce a matrix of transmitted 
disturbances; 

reducing a rank of said matrix of transmitted disturbances to yield a second set 
of basis functions, said second set of basis functions corresponding to composite 
sources, each of said composite sources comprising a linear combination of one or 
more of said original basis functions; 

partitioning a first set of weighting functions into groups, each group 
corresponding to one of said regions, each weighting function corresponding to a 
condition, each of said weighting functions corresponding to an original tester; 

calculating a far-field disturbance received by each of said testers in a first 
group for each of said spherical angles to produce a matrix of received disturbances; 

reducing a rank of said matrix of received disturbances to yield a second set of 
weighting functions, said second set of weighting functions corresponding to 
composite testers, each of said composite testers comprising a linear combination of 
one or more of said original testers; and 

transforming said system of linear equations to use said composite sources and 
said composite testers. 

2. A method of data compression, comprising: 

partitioning a first set of basis functions into groups, each group corresponding 
to a region, each basis function corresponding to an unknown in a system of 
equations, each of said basis functions corresponding to an original source; 

selecting a first plurality of angular directions; 
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calculating a disturbance produced by each of said basis functions in a first 
group for each of said angular directions to produce a matrix of disturbances; 

using said matrix of disturbances to compute a second set of basis functions, 
said second set of basis functions corresponding to composite sources, wherein at least 
one of said composite sources produces a relatively weak disturbance from a portion 
of space around said at least one composite source; 

partitioning a first set of weighting functions into groups, each group 
corresponding one of said regions, each weighting function corresponding to a 
condition, each of said weighting functions corresponding to an original tester; 

calculating a disturbance received by each of said testers in a second plurality 
of angular directions to produce a matrix of received disturbances; 

using said matrix of received disturbances to compute a second set of 
weighting functions, said second set of weighting functions corresponding to 
composite testers, wherein at least one of said composite testers weakly receives 
disturbances from a portion of space relative to said at least one composite tester; and 

transforming at least a portion of said system of equations to use one or more 
of said composite sources and one or more of said composite testers. 

3. The method of Claim 2, wherein said matrix of disturbances is a moment 
method matrix. 

4. The method of Claim 2, wherein said step of using said matrix of disturbances 
to compute a second set of basis functions comprises reducing a rank of said matrix of 
disturbances. 

5. The method of Claim 2, wherein said step of using said matrix of received 
disturbances to compute a second set of weighting functions comprises reducing a rank of said 
matrix of received disturbances. 

6. The method of Claim 2, wherein said disturbance is at least one of an 
electromagnetic field, a heat flux, an electric field, a magnetic field, a vector potential, a 
pressure, a sound wave, a particle flux, a weak nuclear force, a strong nuclear force, and a 
gravity force. 
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7. The method of Claim 2, wherein said first plurality of directions is 
substantially the same as said second plurality of directions. 

8. The method of Claim 2, wherein said regions of space around said at least one 
composite source are far-field regions. 

9. The method of Claim 2, wherein said at least a portion of a region around said 
at least one composite tester is a far-field region. 

10. A method of data compression, comprising: 

calculating one or more composite sources as a linear combination of one or 
more basis functions, wherein at least one of said composite sources produces a 
relatively weak disturbance in a portion of space related to said at least one composite 
source; 

calculating one or more composite testers as a linear combination of one or 
more weighting functions, wherein at least one of said composite testers is affected 
relatively weakly by disturbances propagating from a portion of space around said at 
least one composite tester; and 

transforming at least a portion of a first system of equations based on said basis 
functions and said weighting functions into a second system of equations based on 
said composite sources and said composite testers. 

11. The method of Claim 10, wherein said disturbance is at least one of, an 
electromagnetic field, a heat flux, an electric field, a magnetic field, vector potential, a 
pressure, a sound wave, a particle flux, a weak nuclear force, strong nuclear force, and a 
gravity force. 

12. The method of Claim 10, wherein said basis composite sources comprise 
electric currents. 

13. The method of Claim 10, wherein said composite sources comprise magnetic 
currents. 

14. The method of Claim 10, wherein said composite sources comprise acoustic 
sources. 
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15. The method of Claim 10, wherein said composite sources comprise 
electromagnetic sources. 

16. The method of Claim 10, wherein said composite sources comprise thermal 
sources. 

17. The method of Claim 10, wherein each of said composite sources corresponds 
to a region. 

18. The method of Claim 10, wherein said second system of equations is described 
by a sparse block diagonal matrix. 

19. The method of Claim 18, further comprising the step of reordering said sparse 
block diagonal matrix to shift relatively larger entries in said matrix towards a desired corner 
of said matrix. 

20. The method of Claim 10, further comprising the step of solving said second 
system of equations. 

21. The method of Claim 10, further comprising the step of solving said second 
system of equations to produce a first solution vector, said first solution vector expressed in 
terms of said composite testers. 

22. The method of Claim 21, further comprising the step of transforming said first 
solution vector into a second solution vector, said second solution vector expressed in terms 
of said weighting functions. 

23. A method, comprising: 

calculating at least one composite source, said composite source representing 
energy sources; 

calculating at least one composite tester; and 

transforming at least a portion of a first system of linear equations into a 
second system of linear equations based at least on said at least one composite source 
and said at least one composite tester. 

24. The method of Claim 23, wherein said at least one composite source represents 
a linear combination of one or more energy sources such that said at least one composite 
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source radiates relatively little energy into a portion of angular region disposed about said at 
least one source. 

25. The method of Claim 23, wherein said at least one composite tester is affected 
relatively weakly by energy propagating from a portion of space around said at least one 
composite tester. 

26. The method of Claim 23 , wherein said second system of linear equations is 
represented by a block sparse matrix. 

27. An apparatus comprising: 

means for calculating at least one composite source; 
means for calculating at least one composite tester; and 

means for transforming at least a portion of a first system of equations into a 
second system of equations based at least on said at least one composite source and 
said at least one composite tester. 

28 . A method of data compression, comprising: 

calculating one or more composite sources as a combination of one or more 
basis functions, wherein at least one of said composite sources produces a relatively 
weak product in a portion of space; 

calculating one or more composite testers as a combination of one or more 
weighting functions, wherein at least one of said composite testers interacts relatively 
weakly with said at least one composite tester; and 

transforming at least a portion of a first array of interaction data based on said 
basis functions and said weighting functions into a second array of interaction data 
based on said composite sources and said composite testers. 

29. The method of Claim 28, wherein said disturbance is at least one of, an 
electromagnetic field, a heat flux, an electric field, a magnetic field, vector potential, a 
pressure, a sound wave, a particle flux, a weak nuclear force, strong nuclear force, a gravity 
force, and an image element. 

30. The method of Claim 28, wherein each of said composite sources corresponds 
to a region. 
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31. The method of Claim 28, wherein said second array of interaction data is 
described by a sparse block diagonal matrix. 

32. The method of Claim 28, further comprising the step of using said second array 
of interaction data to compute a first solution vector, said first solution vector expressed in 
terms of said composite testers. 

33. The method of Claim 32, further comprising the step of transforming said first 
solution vector into a second solution vector, said second solution vector expressed in terms 
of said weighting functions. 
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COMPRESSION AND COMPRESSED INVERSION OF INTERACTION DATA 

Abstract of the Disclosure 



A compression technique compresses interaction data. A fast method processes the 
compressed data without the need to first decompress the data. In one embodiment, the 
compression technique is used to compress data in an interaction matrix. The interaction 
matrix (such as a moment method impedance matrix) contains interaction data between 
sources (e.g., basis functions or expansion functions) and testers (e.g., testing functions). The 
sources are collected into groups of sources according to specified criteria. One useful criteria 
is based on grouping sources relatively close to one another. For each group of sources, a 
composite source is calculated. The testers are also collected into groups and composite 
testers are calculated. The use of composite sources and composite testers to compute 
couplings when the source and tester are not close to each other allows the interaction matrix 
to be computed as a sparse matrix with a block format. 
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